Optimization in random field Ising models by quantum annealing 
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We investigate the properties of quantum annealing applied to the random field Ising model in 
one, two and three dimensions. The decay rate of the residual energy, defined as the energy excess 
from the ground state, is find to be eres ~ log(A^Mc)~'' with in the range 2. ..6, depending on 
the strength of the random field. Systems with "large clusters" are harder to optimize as measured 
by ^. Our numerical results suggest that in the ordered phase ( = 2 whereas in the paramagnetic 
phase the annealing procedure can be tuned so that ( ^ 6. 
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I. INTRODUCTION 

It is desirable to have an optimization method which 
can be apphed to as wide range of problems. An example 
method is the classical simulated annealing (SA). During 
the recent years quantum annealing (QA) has gained a 
lot of attention as a promising candidate for a method, 
with a promise of a faster convergence to the optimal 
configuration for a given problem. This has been partly 
motivated by real world realizations, demonstrated in ex- 
periments by Brooke et al. l]. 

Test problems for QA can be found from many 
problem-specific optimization algorithms which find the 
exact ground state in a polynomial time, as for instance 
the random field Ising model or the Ising spin glass in 
two dimensions i^j- A convenient measure for the effec- 
tiveness of the annealing is the residual energy e^es which 
gives the energy difference between the true ground state 
and the configuration that is obtained in the end of the 
annealing process. The most important quantity is the 
annealing time r, during which the temperature (trans- 
verse field in the case of quantum annealing) of the sys- 
tem is reduced to zero. When r is increased the energy of 
the resulting configuration approaches the ground state 
energy. 

For classical simulated annealing it has been predicted 
by Huse and Fisher ^3] that the residual energy e^es de- 
creases with the annealing time r as 



log(T) 



-C 



(1) 



For large time scales they have derived an upper limit 
for the decay of the residual energy in two- level systems: 
C < 2. This result is argued to hold also for random 
field magnets and other disordered spin systems. The 
existence of a phase transition can change the value of 
C- According to Ref. when a random field magnet is 
cooled through the phase boundary from the disordered 
to the ordered phase the annealing slows down to C = 1 . 

In Ref. 4] Santoro et al. have studied quantum anneal- 
ing in the case of a two dimensional (2d) spin glass. They 
have derived a theoretical estimate for the decay rate of 
the residual energy. Santoro et al. argued that the resid- 
ual energy comes from tunnelings at (avoided) Landau- 
Zener (LZ) crossings. The corresponding average resid- 



ual energy resulting from this process was estimated as 
eresir) = drZ(r)Se.(r)exp(r/re(r)) where Z(r) is 
the density of LZ crossings and E^x (F) is the correspond- 
ing average excitation energy. The term exp(T/Tc(F)) 
gives the probability that the system tunnels during the 
annealing to a higher energy eigen-state due to LZ cross- 
ings, and Tc(r) ~ exp(A/£(r)) where ^ is a typical wave 
localization length. Ref. jj] now considers the limit F — > 
which is expected to dominate the annealing behavior. 
With the assumptions ^(F) - F"^ and Z{T)Eex{T) - F'^ 
one gets e^es ^ log(T)^^ with ^ = (1 + uj)/f. With es- 
timates (p — 1/2 and oj = 2 given in Ref. 4] one gets 
^ = 6. The value for tp comes from a quasi-classical 
consideration of a particle's wave length, whereas uj — 2 
can be reasoned for as follows. In the limit F ^ the 
tranverse field can be considered as a perturbation for 
which only the second order correction to the energy has 
a non-zero value. From this follows that Ef.x ~ F^. If 
the density of LZ crossings in the limit F — > is assumed 
to be at most the density of the classical states at F = 
one gets ui = 2. Thus the residual energy decreases log- 
arithmically but now with much larger exponent C ~ 6. 
This implies a considerable speed-up compared to ^ = 2, 
the upper limit of SA. It is useful to emphasize that this 
estimate for is obtained without any problem-specific 
assumptions. 

The numerical results of Santoro et al. (3| show that 
the residual energy of the 2d Edwards-Anderson spin 
glass converges indeed at much faster rate and to lower 
values with QA, compared to SA, though no empirical 
values of ^ were given. In addition, the performance of 
QA has been tested also on the traveling salesman prob- 
lem ^ for which similarly to 2d spin glass it was found 
that QA gives a faster decay of the residual energy than 
SA. However, this is not generally valid for all optimiza- 
tion problems. Battaglia et al. (Qj] have found that in the 
case of the three-satisfiability problem quantum anneal- 
ing is outperformed by simulated annealing. In small sys- 
tems, it is possible to see a power law decay of the residual 
energy, e.g. by solving the time-dependent Schrodinger 
equation adiabatically However, with increasing 

system sizes the energy gap for the Landau- Zcner cross- 
ings decreases 0, which means that the probability of 
staying in the ground state decreases as well leading to 
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the logarithmic behavior discussed above. 

In this paper we study the quantum annealing applied 
to the random field Ising model. The RFIM has the 
following Hamiltonian: 

H = - J^^SiSj -^^hiSi, (2) 

where J > is the coupling constant, Si = ±1 are classi- 
cal spin variables and hi is the random field at site i. One 
of the advantages of RFIM as a test problem is the fact 
that its exact ground state can be found in a polynomial 
time with an efficient graph algorithm from combinatorial 
optimization This allows us to calculate the residual 
energy as the difference between the true ground state 
and the configuration given by any annealing procedure. 

Another feature of RFIM is the fact that the phase dia- 
gram depends on both dimension and the strength of dis- 
order. The second order phase transition of the 2d Ising 
model is destroyed by the random fields, though resid- 
ual ordering persists in finite systems 10]. Though there 
is thus no long range order either in one or two dimen- 
sions, systems of a finite size may have ordered ground 
states when the typical cluster size exceeds the system 
size, which is true also at zero temperature [I^, [HI [l3 ■ 
However, in three dimensions one has a temperature de- 
pendent critical strength of the random field hc{T) be- 
low which the system is ordered '13^. At zero tem- 
perature its value has been calculated numerically as 
hc{T = 0) = 2.27 

We calculate the residual energy as a function of an- 
nealing time measured in Monte Carlo steps in one, two 
and three dimensions with varying strength of disorder. 
Our numerical results suggest that the residual energy 
decays logarithmically as in Eq. However, the value 
of the exponent C, now seems to depend on the nature 
of the ground state. In the ferromagnetic case C « 2, 
whereas in the paramagnetic case C ~ 6, if the cooling 
schedule is tuned well enough. To our knowledge this is 
the first time when empirical values for C, are presented if 
not counting the simple models giving a power law dacay 
of the residual energy 0, |1| . 

The structure of the rest of the paper is the following. 
In section Ull we briefly review the Suzuki- Trotter map- 
ping which transforms a d-dimensional quantum system 
to a d -|- 1-dimensional classical one making the problem 
accessible for conventional Monte Carlo sampling. Sec- 
tion ^H] is devoted to the numerical results starting with 
the results for the classical simulated annealing, which 
serve as a measuring stick when the efiiciency of the 
quantum annealing is discussed later on. For quantum 
annealing we discuss numerical results in one, two and 
three dimensions for different random field strengths. It 
is also studied how the performance of QA is affected 
when the parameters of the d + 1-dimensional equivalent 
classical system are varied. The paper is summarized in 
Section Hfl 



II. NUMERICS 

The quantum version of Eq. jSJ is obtained by replac- 
ing the spin variables Si with Pauli spin operators af. 
Quantum fiuctuations are tuned by changing the strength 
F of a perpendicular field term, arising from the Pauli 
spin operator . 

In the quantum annealing one starts with a large value of 
F so that spins in the z direction are totally uncorrelated. 
By decreasing F gradually towards zero spins fall into the 
ground state configuration provided that T = 0. 

With the Suzuki-Trotter mapping this d- 

dimensional quantum system can be represented by P 
coupled replicas of the classical system, Eq. ^ resulting 
in a d-\- 1 dimensional classical problem with the following 
Hamiltonian: 

= - E I J E + E ^^^"i + E ' 
fc=i \ (jj) I i j 

(4) 

where is the F -dependent coupling constant between 
the replicas: 

PT F 
J_L = Y^^'^^-a^-pf- (5) 

The resulting system has periodic boundary conditions 
in the extra dimension. It is convenient to set periodic 
boundaries also for the original classical system. The 
annealing of the Hamiltonian (Eq. is simulated by a 
standard Monte Carlo sampling of Eq. Q at the effective 
temperature PT with gradually decreasing F. 

The values of random field hi are taken from the fixed 
Gaussian distribution Pcihi) with the parameters (hi) — 
and (hf) — 1. The strength of the random field is tuned 
by varying the ferromagnetic coupling constant J. The 
residual energy e^es is calculated as 

eres{NMc) ^ {Ecl{NMc))k - Egs, (6) 

where {Eci)k is the average energy of all replicas and Eqs 
is the ground state energy. The GS energy and config- 
uration are computed for each sample, as noted in the 
introduction, by using a combinatorial optimization al- 
gorithm. Both {Eci{NMc))k and Eqs from the definition 
in Eq. (|SJ) are normalized per spin. 

III. RESULTS 

A. Classical simulated annealing 

First, we briefly present numerical results for classi- 
cal simulated annealing, which are then later utilized 
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FIG. 1: A doubly logarithmic plot of residual energies in 
classical simulated annealing. L = 32 in 3d, L = 256, in 2d 
and L = 10'* in Id with diflFerent coupling constants J and 
initial temperatures Tq. The straight lines are guides for the 
eye. 



FIG. 2: A doubly logarithmic plot of residual energies in Id 
corresponding to different quantum annealing schedules Eqs. 



in a comparison with quantum annealing. Since the 
differences between various trial annealing schedules in 
the case of SA turned out to be rather small we show 
in Fig. n only the data that corresponds to a linear 
cooling schedule. In one and two dimensions we find 
the expected logarithmic decrease of the residual energy 
feres ~ logfA^Mc)"'') with C, — 2. ..3. According to Ref. 
PI C = 2 is an upper limit, and hence it is expected that 
asymptotically C ^ 2. 

In 3d in the ordered phase, for J > Jc ~ 0.44 [T^ . 
one can observe the slowing effect of the phase transition 
on the cooling efficiency in agreement with the theoret- 
ically predicted C = 1 0- This is clearly visible for the 
case where the annealing is started at low temperature 
(ovals). When the starting temperature is well above Tc 
the system stays in the paramagnetic phase during most 
of the simulation time resulting in ^ = 2. With increas- 
ing number of Monte Carlo steps the system is expected 
to spend more annealing steps in the ferromagnetic phase 
and hence to experience the slowing of the annealing rate 
to C = 1- 



B. Quantum annealing: preliminaries 

A priori it is naturally not clear which is the best way 
to reduce the value of F. We have tested the three fol- 
lowing annealing schedules: 

N 

TuniN) = Fo(l - -— ), (7) 



AT p 

FR(7V)=Fo/(l + i?-— ) , i?=-£-l, (8) 

rLog(iV) = -log{ tanh [ atanh(e-^«) (9) 
- (atanh(e"^°) - atanh(e"^-f )) 1 }. 

The residual energies corresponding to the different 
schedules (Eqs. 0-®) for Id systems of 10"' spins (aver- 
aged over more than 10 samples each) are shown in Fig. 
El With the used set of parameters ( J = 1, PT = 4, P = 
128, Fo — 8,F/ = 10^^) the residual energies seem to 
decay with the same slope for all annealing schedules. 
However, the logarithmic schedule (Eq. 0) gives the 
lowest residual energy and hence we restrict ourselves 
to it throughout the rest of the paper. The choice of 
Fq does not alter the results as far as it is chosen large 
enough in order to ensure that there are no correlations 
in the starting configuration. We flip only one spin at a 
time. We tested also the use of the global flips where one 
attempts to flip all replicas of a given spin. It turned out 
that the single flip strategy is more effective. 

The typical evolution of the quantum annealing in Id 
is illustrated with snap-shots of the spin configurations in 
Fig. 01 The 256 replicas of the original classical system of 
size L — 256 lie in the horizontal direction. Black pixels 
represent spins which match with the orientation in the 
ground state configuration, while white pixels correspond 
to incorrectly aligned spins. Note how the effectively 
classical system is strongly correlated along the Trotter 
direction. 
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FIG. 3: Snapshots of the spin configurations in the process of 
the annealing after (from top to bottom) l/'iNMC, 2/3Nmc 
and Nmc (final) Monte Carlo iterations. Black color de- 
notes correctly aligned spins with respect to the ground state. 
Right: J = 2.0, left: J = 4.0 (L = 256, P = 256, PT = 2, 
To — 40, logarithmic annealing scheme). 



C. Quantum annealing: one dimension 

The Suzuki- Trotter mapping is exact only in the limit 
PT — i- oo. This means that the quantum nature of an- 
nealing should not be seen when the value of PT is too 
small. On the other hand in order to find the ground 
state of the system the temperature T needs to be taken 
to zero. In the ideal case one should have PT — > oo and 
T = 0, simultaneously, which is impossible in practice. 

With P — 1 one simply performs Metropolis dynamics, 
and as P is increased the annealing should become more 
efficient. This is demonstrated in Fig. 0| for Id systems 
of 10* spins. The effective temperature is kept constant 
PT = 4 while the number of the Trotter replicas P is 
varied. The results are averaged over 10 realizations of 
disorder, each. 

When P is increased sufficiently (P > 32) one can 
observe a region where the annealing rate of QA is con- 
siderably higher compared to SA. With increasing P this 
region grows and we expect that this is the asymptotic 
behavior for P — > oo. Since the Suzuki- Trotter mapping 
is exact in this limit, one can assume that this region 



reflects the properties of quantum annealing with a fi- 
nite, non-zero quench rate. For finite values of P the 
system falls out of this quantum annealing regime to the 
classical Metropolis behavior as Nmc is increased. In the 
classical regime the system is fully correlated in the Trot- 
ter direction and the annealing process cannot anymore 
take an advantage of the extra, non-classical dimension. 
This means that in order to maintain a fast annealing 
rate for large values of Nmc larger P values are needed 
as well, as was also observed in Ref. Due to the 

finite simulation temperature the residual energy finally 
saturates to some non-zero value. 

As indicated in Fig. |S1 the decay rate of e^es depends 
on the actual value of the coupling constant J. The case 
with J = and (/i^) = 1 is a problem of L'^ independent 
spins for which the results depend on the used annealing 
schedule ^3 ■ For the logarithmic and rational schedules 
(Eq. lO) we find a polynomial decay of the residual en- 
ergy Cres ~ {Nmc)~'' with « 2 whcrcas in the case 
of the linear annealing schedule ( ^ 1. For J > the 
quantum annealing goes over to the logarithmic regime 
Eres ~ log(A'^A/c)~^- The numerical values of C roughly 
agrees with the estimate Cmax — 6 given by Santoro et 
al. 3. As J grows, and hence the cluster sizes of the 
ground state, the annealing efficiency seems to diminish. 
When no random field is applied the dynamics of the 
quantum annealing in the limit F — > for P » L cor- 
responds to a case of a strongly anisotropic Ising model. 
In Ref. Ferreira et al. have studied the two dimen- 
sional anisotropic Ising system with » ksT » Jy. 
They found that in the large time limit the width of the 
interfaces perpendicular to x-direction saturates to some 
finite value. We have verified numerically that in same 
limit the quantum annealing ends up in a situation where 
the residual energy comes from rough, fluctuating bound- 
aries, positioned perpendicular to the Trotter direction. 

The scale of the energy barriers is determined by the 
value of J. Hence, when J is increased one needs larger 
PT to overcome the barriers. The data in Fig. shows 
how the annealing rate grows as PT is increased. The 
value of PT has a two-fold effect. As indicated in Fig.|H| 
large values of PT are desirable in order to minimize the 
error in the Suzuki- Trotter mapping and hence to be able 
to observe the true quantum annealing behavior. On the 
other hand the PT determines the growth rate of J± (see 
Eq. O . As PT increases one needs larger a Nmc in order 
to keep the same annealing rate and hence to reach an 
equally low value of e^es ■ 



D. Quantum annealing: 2d and 3d 

In two dimensions we consider systems with sizes only 
up to 128 X 128 in order to be able to consider P up to 
1024. For low values of J, where the size of the clusters is 
well below the system size, the QA performs similarly to 
the one dimensional case giving C w 6 for J = 0.33 (see 
Fig. [71). With J also the cluster sizes of the ground state 
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FIG. 4: A doubly logarithmic plot of residual energies in Id 
with J = 1,L = 10*, PT = 4, P = 8. ..1024. For comparison 
we also show the data corresponding to the classical simulated 
annealing (SA). In order to compare the used Monte Carlo 
time, in the case of QA Nmc has to be multiplied with the 
number of replicas P. 



FIG. 6: A doubly logarithmic plot of residual energies in Id 
with L = 10*, J = 1. The annealing rate increases with PT 
at the cost of growing amplitude (prefactor). 
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5: A doubly logarithmic plot of Id residual energies, 



FIG 

L = 10*, PT = 4. The J = case is a problem of L 
independent spins for which eres ~ {Nmc)~'' with ~ 2. For 
J > Crea seems to decay logarithmically. With increasing J 
the RFIM becomes more difficult to anneal. 



grow reaching the system size (128 x 128) approximately 
at J = 1.5. From Fig. □ one can see that for J ^ 1 
one has C = 2. Whereas in the case of one dimensional 
systems with extremely weak disorder the value of C could 
be raised by increasing the effective temperature PT this 
seems not to work in two dimensions. This suggests that 
there is a fundamental difference in the performance of 




log( N 



MC 



FIG. 7; A doubly logarithmic plot of residual energies in 
2d, PT — 12, L — 128, and one data set corresponding to 
L — 16. The straight lines are guides for the eye. For J > 1 
the residual energy decays with « 2. With growing J no 
further decrease of is observed. For comparison we also 
show the results for classical simulated annealing (SA). 



QA depending whether the system has a disordered or 
ordered ground state. This conclusion is supported by 
the observation that with J = 1 for a larger system (L = 
128, squares in Fig. one gets lower residual energies 
compared to a smaller system {L = 16, triangles in Fig. 
\7\) that already has an ferromagnetically ordered ground 
state. 

As it was already evident in the 2d case also the results 
for 3d RFIM show that QA is sensible to the ordering of 
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FIG. 8: A doubly logarithmic plot of residual energies in 
3d, L = 24, J = 0.33. With increasing PT the effective ( is 
growing. For comparison we also show the results for classical 
simulated annealing (SA). 
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FIG. 9: A doubly logarithmic plot of residual energies in 
3d, L — 24, J = 0.66. The decay rate of Cres does not vary 
with PT, ^ ~ 2. For comparison we also show the results for 
classical simulated annealing (SA). 



the underlying system. For low values of J (J < Jc ~ 
0.44 where the system is in the paramagnetic phase 
also at T = we find the same behavior as in one and two 
dimensions (Fig. |SJ): QA is faster than SA and C can be 
tuned towards 6 by increasing PT. Fig.inishows the data 
corresponding to J = 0.66. For the range of parameters 
that have been used we find that C « 2, as in 2d, for a 
system with an ordered ground state. 

IV. SUMMARY 



We have studied numerically quantum annealing in the 
random field Ising model and compared our results with 
classical simulated annealing. When the system is in the 
paramagnetic phase we find that asymptotically QA pro- 
vides a better decay rate of the residual energy with C up 
to 6 in agreement with the Landau-Zener -picture based 
scaling argument presented by Santoro et al. Q. 

We expect that the asymptotic performance of QA 
in one and two dimensions does not change by varying 
the coupling constant J or the magnetic field strength 
h. With growing cluster sizes in the GS one needs in- 
creasingly larger values of PT and P for which the fast 
annealing rate with C « 6 could be observed. The re- 
quirement of large values of P in the case of weak random 
field makes QA from the practical point of view slower 
compared to SA. Note the starting point however: that 
the RFIM GS can be found effectively with combinatorial 
optimization. 

In 3d we have presented evidence that the performance 
of QA depends on whether the system is in the param- 
agnetic or ordered phase. Thus, the situation is actually 
analogous to the behavior of the SA. In the paramag- 
netic phase we find the similar behavior as in one and 
two dimensions. In the ordered phase we observe a much 
slower decay of e^es with ^ « 2, so that in fact QA is 
slower than SA with a starting temperature Tq > Tc- 

We conclude with the general observation that the bet- 
ter efficiency of QA is most clear when the ground state 
consists of small clusters, i.e. the correlation length of 
the ground state is short compared to the number of the 
Trotter replicas. Such benefits vanish with an increasing 
correlation length, of the ground state configuration. 
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